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TRANSIENT  TRANSPORT  IN  NANOSTRUCTURE 
DEVICES  VIA  THE  QUANTUM  LIOUVILLE 
EQUATION  IN  THE  COORDINATE  REPRESENTATION 


ABSTRACT 

Through  the  use  of  numerical  methods  the  quantum  Liouville  equation  in  the  coordinate 
representation  has  been  implemented  for  the  study  of  semiconductor  devices  with  nanoscale  feature 
sizes,  with  an  emphasis  on  dissipation.  We  have  successfully  applied  this  algorithm  to  the 
calculation  of  current  in  single  and  multiple  barrier  diodes,  and  Scientific  Research  Associates,  Inc, 
is  presently  one  of  the  few  research  organizations  numerically  obtaining  the  quantum  distribution 
functions  for  both  electrons  and  holes.  This  document  summarizes  work  performed  under  US 
Navy,  Office  of  Naval  Research  Contract:  N00014-93-C-0044. 
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1.  INTRODUCTION 


In  the  past  decade  we  have  seen  demonstrations  that  band  structure  engineered  devices  with 
nanoscale  dimensions  have  the  potential  of  being  the  basis  for  a  dramatic  shift  in  semiconductor 
technology.  One  significant  illustration  of  this  has  been  the  development  of  multiple  value  logic  and 
memory  circuits  based  on  utilization  of  the  negative  differential  resistance  properties  of  either 
isolated  resonant  tunneling  structures  or  embedded  RTDs  (as  in  the  case  of  resonant  tunneling 
bipolar  transistors). 

Generally,  the  approach  to  studying  these  band  structured  engineered  device  falls  into  two 
categories.  The  first  category  is  concerned  with  studying  fundamental  device  behavior  and 
predicting  device  performance  through  computation  of  the  quantum  distribution  function.  This 
approach  emphasizes  the  details  of  transient  behavior,  the  numbers  of  particles  involved  in  device 
operation,  the  temporal  duration  under  which  the  effective  mass  approximation  is  valid,  the 
significance  of  the  Fermi-golden  rule,  and  other  short  time  phenomena.  In  this  approach,  when 
scattering  is  present,  or  when  time  dependent  fields  are  present  and  treated  as  perturbations,  it  is 
supposed  that  the  perturbation  does  not  modify  the  states  of  an  unperturbed  system,  rather  the 
perturbed  system  instead  of  remaining  permanently  in  one  of  the  unperturbed  states  is  assumed  to 
be  continually  changing  from  one  to  another,  i.e.,  undergoing  transitions  from  one  state  to  another 
state.  This  approach  is  at  the  heart  of  those  calculations  employing  the  density  matrix,  those 
employing  the  Wigner  distribution  function,  and  those  employing  Green's  function  techniques.  The 
second  category,  which  is  presently,  the  most  popular,  obtains  solutions  to  the  time  independent 
Schrodinger  equation  under  nonequilibrium  conditions,  from  which  the  time  independent  wave 
functions  are  obtained.  This  second  procedure  does  not  yield  the  quantum  distribution  fiinction,  and 
cannot  reliably  predict  device  performance.  The  approach  taken  at  Scientific  Research  Associates, 
Inc.,  (SRA),  has  been  to  obtain  the  quantum  distribution  function  as  solutions  to  the  density  matrix 
in  the  coordinate  representation. 

This  document  summarizes  OMR  Contract  N00014-93-C-0044  sponsored  studies  of  the 
density  matrix  at  SRA. 

2.  SUMMARY  OF  STUDIES  UNDER  ONR  CONTRACT  N00014-93-C-0044 

The  density  matrix  studies  at  SRA  have  been  developed  along  very  general  lines.  The  SRA 
code  is  capable  of  treating: 

•  large  numbers  of  particles  present  in  the  transport  process, 

•  the  effects  of  the  Pauli  exclusion  principle  through  boundary  conditions, 

•  electron  and  hole  transport, 

•  and  transients. 

The  algorithm  implemented  has  deliberately  ignored  regional  approximations,  where  the 
quantum  mechanical  region  is  connected  to  the  classical  region  via  boundary  conditions.  Instead 
SRA  has  obtained  the  quantum  distribution  function  for  the  entire  structure  being  studied. 
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2a.  THE  ROLE  OF  DISSIPATION 


The  stnicture  that  we  typically  dealt  with  incorporated  the  specific  quantum  mechanical 
barrier/well  region  and  the  cladding  regions  surrounding  the  barriers.  From  the  outset  of  the  study 
we  incorporated  dissipation.  This  introduces  its  own  set  of  difficulties.  First,  as  was  learned  fi’om 
classical  device  studies,  incorporation  of  dissipation  is  essential  to  any  development  of  a  consistent 
set  of  device  physics.  Dissipation  is  necessary,  if  current  is  to  flow  within  the  device  under  steady 
state  conditions.  Dissipation  is  not  required  for  transient  current  contributions.  The  question 
we  grappled  with  was  what  was  the  best  means  to  incorporate  dissipation.  Clearly,  from  an 
analytical  point  of  view  dissipation  is  best  handled  by  incorporating  it  in  the  governing  equations. 

The  simplest  form  of  dissipation  is  that  associated  with  the  relaxation  approximation.  While 
questions  arise  as  to  whether  the  relaxation  time  should  be  energy  dependent,  and/or  momentum 
dependent,  and/or  position  dependent,  this  is  the  most  frequently  used  approximation.  Indeed 
many  of  us,  who  have  been  involved  with  classical  device  transport  calculations,  routinely 
incorporated  the  relaxation  time  approximation  in  terms  that  arose  from  the  linearized  Boltzm^ 
equation.  The  linearized  Boltzmann  equation  is  often  used  as  the  basis  for  the  drift  and  difiusion 
equations.  But  the  use  of  the  relaxation  time  approximation  does  not  deal  with  one  of  the  naost 
basic  of  device  issues,  the  representation  of  the  regions  between  the  quantum  mechanical  barriers 
and  wells,  and  the  physical  contacts.  It  is  this  region  where  the  carriers  undergo  sufficient 
scattering  to  relax  to  local  steady  state  values,  and  result  in  the  concept  of  a  quasi-Fermi  level.  This 
problem  was  dealt  with  as  discussed  below. 

The  presence  of  a  relaxation  time  in  the  equation  of  motion  of  the  density  matrix  results  in  a 
governing  equation  that  is  complex,  i.e.,  dissipation  introduces  an  imaginary  contribution  to  the 
steady  state  Liouville  equation.  On  the  other  hand  the  introduction  of  a  quasi-Fermi  energy  to  the 
quantum  Liouville  equation,  without  regard  to  a  local  relaxation  time,  results  in  a  net  change  in  the 
local  value  of  the  potential  energy,  without  an  imaginary  contribution  to  the  governing  equation. 
If  these  contributions  are  treated  separately,  then  when  using  the  quasi-Fermi  energy 
approximation,  current  is  introduced  through  boundary  conditions,  through,  e.g.  the  density  matrix 
equivalent  of  a  ^placed  Fermi  distribution. 

At  SRA  a  wide  variety  of  dissipation  models  were  introduced  into  the  simulations.  We 
implemented  the  relaxation  time  approximation  but  did  not  emphasize  its  contribution.  Instead  we 
incorporated  the  concept  of  the  quasi-Fermi  energy,  whose  spatial  distribution  was  determined  by 
the  net  current,  the  local  carrier  density  and  the  local  relaxation  time.  Current  was  introduced 
through  the  presence  of  boundary  conditions  -  in  the  case  of  solutions  to  the  steady  state  quantum 
Liouville  equation.  We  also  introduced  dissipation  through  Fokker-Planck  contributions.  In  the 
case  of  the  density  matrix  this  contribution  takes  the  form: 


1 

/7 

T 
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Whereas  in  the  case  of  Wigner  transport  its  form  is; 
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The  Fokker-Plank  terms  introduce  an  imaginary  contribution  to  the  Liouville  equation  in  the 
coordinate  representation,  and  under  conditions  where  the  local  density  is  constant  the  dissipation 
forces  current  to  flow  with  parameters  varying  according  to  Ohms  Law.  (The  normalizing  term  in 
the  above  equation  is  the  thermal  deBroglie  wavelength.) 

One  of  the  important  aspects  of  this  study  was  an  attempt  to  combine  the  concept  of  the 
quasi-Fermi  level  and  dissipation  as  expressed  through  the  Fokker-Plank  terms.  We  succeeded  in 
carrying  this  out  for  classical  systems,  and  performed  several  successful  numerical  simulations. 
Work  was  begun  on  systems  containing  quantum  mechanical  elements  and  simulations  were 
performed  for  modest  quantum  mechanical  barrier.  By  the  end  of  this  contract  period  robust 
quantum  mechanical  simulations  with  both  Fokker  Plank  and  quasi-Fermi  energy  contributions 
were  not  completed. 

The  status  of  the  SRA  density  matrix  algorithm  at  the  end  of  this  study  is  represented  by  the 
Pi’s  lecture  notes  distributed  at  the  NATO  ASI  on  Quantum  Transport  in  Ultrasmall  Devices  (II 
Ciocco,  Italy,  17-30  July  1994),  a  copy  of  which  is  included  with  this  report.  Supplementing  the 
lecture  notes  were  a  series  of  lectures,  in  which  a  variety  of  examples  illustrating  the  effects  of  delta 
doping,  the  significance  of  the  quantum  potential,  and  initial  transient  calculations.  Additionally, 
the  study  of  dissipation,  as  represented  in  the  SRA  algorithm  was  discussed. 

We  note  that  additional  discussions  of  the  density  matrix  and  the  quantum  potential  are 
contained  in  the  following  papers,  were  published: 

1.  Modeling  of  Quantum  Transport  in  Semiconductor  Devices,  (with  D.  K.  Ferry) 
Published  by  Academic  Press,  Inc. 

2.  Density  Matrix  Simulations  of  Semiconductor  Devices,  Published  by  Plenum  Press. 

3.  Transport  in  Two-Dimensional  Quantum  Well  HEMTs  (with  J.  P.  Kreskovslq') 
Proceedings  of  Third  International  Conference  on  Computational  Electronics. 

4.  Resonant  Tunneling  Calculations  via  the  Density  Matrix  in  the  Coordinate 
Representation  (with  T.  R.  Govindan)  Proceedings  of  Third  International  Conference 
on  Computational  Electronics. 

5.  Electron  Transport  Using  the  Quantum  Corrected  Hydrodynamic  Equations  (with  J.  P. 
Kreskovsky)  Invited  paper  in  VLSI. 

2b.  NON-SELF  CONSISTENT  TRANSIENT  STUDIES 

The  task  associated  with  the  transient  calculations  was  to  develop  a  transient  accurate 
algorithm  within  the  framework  of  the  quantum  Liouville  equation  in  the  coordinate 
representations.  As  in  the  case  of  the  steady  state  calculations,  the  transients  were  performed  along 
‘characteristic  directions’.  Thus,  much  of  the  effort  here  was  in  algorithm  development, 
implementation  and  testing. 


As  in  most  numerical  simulations  involving  transient  there  are  two  time  intervals  to  deal 
vAth.  The  first  time  interval  is  the  physical  time  interval,  the  second  time  interval  involves  the 
number  of  iterations  required  to  get  to  a  converged  solution  within  the  prespecified  time  interval. 
During  the  initial  stages  of  the  transient  algorithm  development  we  focused  exclusively  on  the 
quantum  Liouville  equation,  ignoring  Poisson’s  equation.  We  examined  both  closed  and  open 
systems,  and  calculations  were  undertaken  that  demonstrated  the  efficacy  of  the  algorithm.  We 
illustrate  the  non-self  consistent  calculations  below. 

The  problem  we  considered  was  for  a  closed  system.  The  initial  state  of  the  system 
consisted  of  an  initial  distribution  of  charge  at  the  left-hand  side  of  the  device.  While  we  could 
have  specified  an  initial  charge  distribution  in  the  form  of  a  Gaussian  —  we  didn  t.  Rather  we 
computed  a  charge  distribution  fi’om  a  potential  well  configuration.  The  position  of  the  potential 
well  used  for  the  initial  charge  distribution  is  displayed  in  figure  1,  where  it  is  represented  by  the 
t  <  0  data.  The  position  of  the  potential  well  used  in  the  transient  calculation  is  represented  by  the 
/ >  0  data.  This  t>0  potential  well  represents  an  attractive  force  on  the  carriers.  The  steady  state 
charge  distribution  arising  from  this  potential  well  is  shown  in  figure  2,  and  is  represented  by  the 
t  =  0  data  points.  We  reiterate:  the  initial  potential  distribution  has  no  impact  on  the 
computations  for  t  >0 .  The  initial  potential  distribution  is  not  present  at  times  r  >  0 .  It  is 
introduced  solely  to  provide  the  initial  charge  distribution. 

The  initial  charge  distribution  consisted  of  a  buUdup  of  charge  at  the  left  side  of  the 
structure.  Subsequent  to  this,  see  figure  2,  is  a  buildup  of  charge  over  the  right-side  quantum  well. 
The  structure  of  the  space  charge  distribution  reveals  that  there  are  strong  nonuniformities  in  the 
charge  distribution  situated  around  the  initial  distribution.  We  know  of  course,  firom  the  discussion 
of  the  Bohm  quantum  potential  that  there  are  quantum  mechanical  forces,  whose  strength  is 
determined  by  the  shape  of  the  carrier  distribution.  (We  note  that  a  similar  interpretation  can  be 
extracted  fi'om  calculations  where  we  watch  the  evolution  of  a  distribution  of  weighted  plane  wave 
states.)  The  implication  of  this  result  is  that  the  relaxation  to  steady  state  is  determined,  in  part, 
by  memory  of  the  initial  charge  distribution.  The  physical  time  sequences  are  indicated  in  the 
figure  2,  and  we  begin  to  see  that  the  transition  to  steady  state  appears  to  involve  a  damped 
oscillation  in  the  charge  distribution,  about  its  initial  and  final  distribution. 

There  are  two  aspects  of  this  calculation  that  are  noteworthy.  First  is  the  demonstration 
that  there  are  memory  effects  and  that  these  damp  out.  Second,  is  the  length  of  the  structure  used 
to  perform  the  calculations.  The  device  structure  is  200  nm  long.  Most  quantum  device 
simulations  are  for  structures  significantly  smaller  than  200  nm.  These  cannot  separate  the 
quantum  mechanical  region  from  the  influence  of  the  boundary. 
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Figure  1.  Initial  and  final  states  associated  with  the  non-self  consistent  time  dependent 
calculation. 


Figure  2.  Time  dependent  distribution  of  charge  associated  with  the  placement  of  the 
quantum  well  to  the  right  hand  side  of  the  structure. 
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2c.  SELF-CONSISTENT  TRANSIENT  STUDIES 

At  the  end  of  the  Contract  period  we  succeeded  in  demonstrating  that  the  transient  self- 
consistent  calculations  could  be  performed.  The  calculation  was  performed  for  a  single  barrier 
structure  in  which  a  change  in  bias  occurred.  Unlike  the  non-self  consistent  studies,  these  dealt 
with  open  boundary  conditions  in  which  the  time  dependent  potential,  density  and  current  were 
calculated 


2d.  ANALYTICAL  STUDIES 

The  time  dependent  density  matrix  calculations  include  terms  that  are  consistent  with  the 
Fermi  golden  rule.  As  such,  while  they  contribute  to  femtosecond  scale  transients  they  do  not 
incorporate  the  necessary  physics  associated  with  short  time  transient.  While  subsequent  work  at 
SRA  in  other  programs  have  attempted  to  generalize  dissipation  and  transients  on  a  short  time 
scale,  the  work  here  involved  teasing  the  density  matrix  algorithm  to  separate  out  explicit  short 
time  and  long  time  scale  events.  This  was  accomplished  by  submitting  the  density  matrix  to  second 
order  perturbation  theory. 

The  second  order  contribution  to  the  time  dependent  equation  of  motion  in  the  coordinate 
representation  is: 


(1) 
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In  the  above  expression  the  term  |M(k,k-q)|2  represents  the  square  of  the  matrix  element  associated 
with  the  perturbation,  and  the  second  contribution  in  the  large  curly  brackets  represents  differences 
in  the  zeroth  order  density  matrix  between  different  quantum  states.  As  a  sanity  check  we  note  that 
in  the  long  time  limit  the  above  expression  yields  standard  perturbation  theory  with  the  Fermi 
golden  rule. 

More  generally  after  performing,  formally,  the  necessary  time  integrations,  the  density 
matrix  to  second  order  in  perturbation  theory  provides  the  following: 
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The  above  equation  is  not  very  transparent,  but  if  we  perform  the  indicated  time  integration, 
the  first  two  contributions  provide  an  indication  of  the  early  time  transients,  which  are  shown  in 
equation  (3).  Equation  (3)  indicates  that  an  expansion  of  the  exponential  contribution  will  yield  and 

early  time  contribution  that  is  of  order  .  The  significance  of  this  result  is  that  it 

indicates  that  it  may  be  possible  to  break  up  the  scattering  contribution  into  explicit  short 
time  and  long  time  contributions. 
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3.  SUMMARY  AND  COMMENTS  ON  RECENT  STUDIES 

The  discussion  of  the  previous  sections  summarizes  the  work  completed  under  Contract 
N00014-93-C-0044.  What  are  some  of  the  issues  that  are  germane  to  this  study  that  have 
transpired  since  this  study  was  completed? 

First,  the  spate  of  device  simulations  utilizing  either  the  density  matrix  in  the  coordinate 
representation  or  the  Wigner  function  analysis  have  provided  us  with  demonstrations  of  what  these 
simulations  can  do.  But  these  simulations  have  not  been  used  for  detailed  device  design,  and 
this  is  a  source  of  disappointment.  (The  closest  attempt  at  a  combined  theory  and  experiment 
study  has  been  that  of  the  NEMO  team,  but  their  study  has  not  produced  a  transient  analysis  and 
there  are  other  issues  to  be  dealt  with  as  we  discuss  below.)  Indeed,  with  the  exception  of  the 
density  matrix  in  the  coordinate  representation,  undertaken  by  workers  at  SRA,  no  other  quantum 
device  simulation  has  reproduced  Ohms  law! 

Analyzing  the  different  approaches  of  these  device  simulations,  in  particular  the  Wigner 
simulations,  indicates  both  physics  and  numerical  based  deficiencies.  The  latter  is  associated  with 
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the  number  of  grid  points  used  for  the  study,  the  size  of  the  device  simulated,  and  the  algorithm 
used  to  convert  the  governing  partial  differential  equation  to  the  required  difference  equation. 

What  are  some  of  the  physics  based  issues  particularly  as  they  affect  device  physics  and 
dissipation?  As  discussed  above  the  ones  of  interest  include  the  role  of  dissipation.  We  know  from 
classical  physics  that  a  solution  to  the  collisionless  Boltzmann  equation  does  not  yield  a  finite 
steady  dc  current,  although  there  can  be  a  transient  current.  Similarly  the  collisionless  Wigner 
equation  does  not  yield  a  steady  dc  current,  only  a  transient  current  without  the  incorporation  of 
dissipation  in  the  governing  equation. 

We  point  out  that  advances  in  classical  device  simulation,  which  include  Gunn  devices, 
IMPATTs,  HBTs,  FETs,  do  not  treat  current  as  a  boundary  condition — it  is  an  integral  part  of  the 
simulation’  Thus  when  we  see  published  papers  of  simulations  with  and  without  scattering  we  need 
to  ask  what  they  signify  and  how  they  can  be  performed  properly.  The  suggestion  for  a  suitable 
approach  is  as  follows:  The  device  simulated  should  consist  conceptually,  of  three  distinct  regions, 
two  boundary  regions  probably  each  50nm  to  lOOnm  long  surrounding  the  third  region  containmg 
the  quantum  structure.  The  boundary  regions  should  be  dissipative,  as  they  would  be  in  an  actual 
structure.  The  quantum  mechanical  region  can  initially  be  modeled  ignoring  scattering  assuming 
that  all  transport  is  ballistic.  Then  for  comparison,  scattering  is  introduced  into  the  quantum 
mechanical  region. 

The  heart  of  the  matter  is  this:  after  the  quantum  mechanical  nature  of  device  transport 
is  established  the  role  of  dissipation  must  be  explicitly  incorporated  into  the  analysis. 
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